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PREFACE 


The revitrf was conducted during a four week visit by Dr. Ken Morgan, 
University College of Swansea, University of Wales, U.K. , to the NASA/ 
Langley Research Center in early summer of 1982. This assessment is based 
primarily upon his and the principal investigator's knowledge of current 
finite element literature. In addition, several important technical discus- 
sions were held with NASA researchers during Professor Morgan's visit, and 
the authors would like to express their appreciation for these helpful 
discussions. 
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FINITE ELEMENT ANALYSIS OF AERODYNAMIC HEATING IN THREE- 
DIMENSIONAL VISCOUS HIGH SPEED COMPRESSIBLE FLOW: 

AN ASSESSMENT 

By 

Ken Morgan^ and Earl A. Thornton^ 

SUMMARY 

The current capability of the finite element method for solving 
problems of viscous flow is reviewed. Much work has been directed to the 
simulation of incompressible flows and the relevant features are described. 
The methods available for, and the problems associated with, the finite 
element solution of high speed viscous compressible flows are analyzed. 

A plan for developing finite element research in this area with experimental 
support is presented. 


INTRODUCTION 


Thermal stresses and deformations induced by aerodynamic heating on 
advanced space transportation vehicles are important concerns in structural 
design. Nonuniform heating may have a significant effect on the performance 
of the structures, and effective analytical ...^thods for predicting the 
structural response are required. For the past few years, the principal 
investigator has been working closely with the NASA/Langley Research Center 
in developing new finite element methodology for thermal-structural analy- 
sis. The finite element method has excellent capabilities for stress analy- 
sis of complex structures, but its capabilities for heat transfer and flow 
analysis are less well-developed. Some recent progress has been made in 
development of finite element methodology for heat transfer, although much 
remains to be done before the method is developed to its full potential. 

1 Lecturer, Department of Civil Engineering, University College of Swansea, 
Wales, U.K. , SA28PP. 

^Associate Professor, Department of Mechanical Engineering and Mechanics, 

Old Dominion University, Norfolk, Virginia 23508. 



For aerodynamic flow analysis the method is still relatively undeveloped in 
comparison to the mature state of the finite difference nwthod in Computa- 
tional Fluid Dynamics (CFD). The finite eleir?nt method is not competitive 
with popular CFD finite difference techniques such as MacCormack's algorithm 
(refs. 1-3). Yet the finite element method may offer significant 
computational advantages in comparison to existing CFD methods and deserves 
further investigation. This report presents a review of recent finite 
element progress in viscous flow analysis with recommendations for the 
important research required to extend the method to determine the 
aerodynamic heating on space transportation vehicles in three-dimensional 
viscous compressible flows. 

VISCOUS INCOMPRESSIBLE FLOW SIMULATION BY 
THE FINITE ELEMENT METHOD 

Two-Dimensional Studies 

Consideration will be restricted to those analyses which have utilized 
the primitive variable formulation in which the basic unknowns are the velo- 
city components and the pressure. The powerful stream-function/vorticity 
methods will not be considered as this approach is not directly applicable 
to the analysis of three-dimensional problems which is our ultimate aim. 

The initial work of Taylor and Hood (ref. 4) used the steady state 
equations with an 8-noded element for the velocity field. Numerical experi- 
mentation soon indicated serious problems when the pressure was approximated 
in a similar manner and they were led to the use of 4-noded bilinear ele- 
ments for the pressure. The presence of the convective terms in the momen- 
tum equations leads to a non-syrametric non-linear matrix system and Hood 
(ref. 5) developed a non-symmetric frontal solver to obtain the solution 
directly at each stage of a standard iteration process. This direct 
approach to the solution has been adopted in the majority of the subsequent 
work in this area. Other element types have been implemented (ref. 6), but 
generally the restriction to a mixed interpolation for velocity and pressure 
applies. 

If the pressure is not of direct interest, it can be effectively 
removed from the analysis (with a consequent reduction in computer core and 
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cose) by adopting the penalty function foriaulation (ref. 7). The pressure, 
p, is approximated by 

p«-Xdivu_ (1) 

where u is the velocity vector and X is a large positive nunA>er termed 
the penalty parameter. The standard incompressibility constraint is then 
replaced by equation (1) and the pressure can be removed from the equation 
system. This approach works well in practice and the convergence of the 
solution to the incompressible solution X ♦ • is now the subject of 
theoretical study (ref. 8). 

In the solution of transient problems, the usual approach adopted has 
been a straightforward implicit algorithm with variants of the Newton- 
Raphson method to speed convergence (ref. 6). A notable exception is the 
work of Donea et al. (ref. 9) who produced a finite element explicit frac- 
tional step method for time-dependent problems that was based on an algo- 
rithm originally developed in the context of finite differences by Qtorin 
(ref. 10). At each time step an auxiliary velocity field is first computed, 
which accounts for all contributions to the momentun equations, except those 
arising from the pressure. The pressure field is then ob^ained by solving 
the Galerkin equivalent of a Laplace type equation, and the time step is 
comple._ed by adding pressure contributions to the auxiliary velocities to 
ensure maintenance of incompressibility. 

Over the past few years, the literature in this area has been full of 
discussions about the need for upwinding. This is a stabilization technique 
which has been used by many investigators in problems where the convective 
terms dominate (ref. 11). The basic idea is to compensate for numerical 
instability due to convection by adding a certain amount of numerical dis- 
sipation, Figure 1. It is now generally agreed that this added dissipation 
should be anisotropic with a component only along the velocity vector (ref. 
12). Some researchers, notably Gresho et al, (ref. 13), have argued against 
the use of upwinding techniques, preferring instead mesh refinement in 
"problem" areas of the flow field. The streamline upwinding method does 
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involve a free parameter whose value may be exactly determined for steady 
linear problems in one dimension but generally is to be selected so as to 
maximize the accuracy. In a recent paper, Donea (ref. 14) claims to have 
overcome some of the problems associated with the analysis of transient 
convective transport processes, and he identifies the source of the trouble 
for hyperbolic problems to be the accuracy of the difference approximation 
to the time derivative term. He shows that a specific cure has to be 
devised for each time-'integration method and the results obtained for simple 
problems do appear to be an improvement over those produced by conventional 
approaches (eg. Figure 2). A detailed assessment of the impact of the tech- 
nique to more complicated problems is required, however, before specific 
conclusions may be drawn. 

Turbulence effects are also being included in finite element models, 
and steady turbulent flow has been analyzed successfully using both 
algebraic (ref. 15) and more sophisticated (ref. 16) turbulent viscosity 
models. 


Three-Dimensional Studies 

This is a new area for the finite element method, and it is interesting 
to notice how it is developing. Tne major contribution so far has been made 
by the group of Gresho at Lawrence Livermore Laboratory (refs, 17-19) which 
is now reporting the solution of problems involving 6400 elements with 
approximately 45,000 equations. These would appear to be the largest finite 
element flow computations ever undertaken. To perform these computations 
efficiently on a CRAY computer, Gresho has left the implicit formulation, 
which he favoured in two-dimensions, and has adopted an explicit scheme of 
the Chorin/Donea type and with the simplest elements. He is also investi- 
gating the effects of using a highly vectorizable one point quadrature to 
evaluate the Galerkin integrals instead of the standard 2’<2^2 Gauss point 
distribution which is not so vectorizable. Initial results showed no signi- 
ficant variation in accuracy between the two approaches, but a substantial 
saving in computer cost is achieved from 7 sec CPU and 13 sec I/O per time 
step to 0.3 sec CPU and 1.3 sec I/O per time step. 'Ihe work of this group 
has indicated, as might be expected, that major problems with computer 





storage and execution time can be encountered in the simulation of coiqplex 
three-dimensional flows, but they have also shown ways in which these 
problems may be overcome and high quality results can be produced. 

Reddy (ref. 20) is also attempting three-dimensional flow calculations 
and is employing a three-dimensional form of the penalty method. He reports 
calculations for natural convection in a cubical box, but admits that his 
computational facilities are such that he is unable to refine the mesh to 
investigate the accuracy of the solution. 

VISCOUS COMPRESSIBLE FLOW SIMULATION BY 
THE FINITE ELEMENT METHOD 


The most significant contribution in this area is the work of Baker 
(refs. 21-23) who has developed an approach for analyzing three-dimensional 
viscous compressible flows. He uses a "dissipative" finite element model 
which, when solving 

L(*) - 0 in n (2) 


replaces the classical weighted lesidual statement 


/ W.L(^) dO - 0 (3) 

£1 ^ 


by 


/ W.L(4») dfl + 6 / W. V(l(*)] d£l - 0 (4) 

£1 ^ £1 ^ 

By performing the usual stability analysis on a linearized one-dimensional 
equation. Baker is able to show that with appropriate choice of 0 the 
resulting linear finite element algorithm is more accurate than an equiva- 
lent finite difference form in which the dissipation takes the form of a 
standard artificial viscosity. The region of interest is initially mapped 
into the unit cube, and the Jacobian of the transformation is interpolated 
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over each element in terms of element nodal values. The quantities J~ ^ 

and det J are stored for each node of each element, and the need for 
-e 

numerical integration during the running of the program is then removed by 
performing a number of calculations beforehand, and storing the relevant 
"hypermatrix" information. The final non-linear matrix equation is solved 
by a direct Newton-Raphson procedure, invoking the tensor matrix product 
method to significantly reduce the core and CPU requirement. Ihis approach 
corresponds to the method of approximate factorization trhich is widely used 
with the finite difference method and consists of replacing the left-hand 
side matrix A in the Newton-Raphson procedure by 

(5) 

where, for example, is constructed in the same raa''.aer as A but con- 

sidering only one diroensioncl interpolation and differentiation in the 
transformed direction S. 

The model has been "tuned" by examining its performance on the one- 
dimensional shock tube problem, but only computations with the parabolized 
form of the equations appear to have been reported in three dimensions. 

The work of Cooke (ref. 24) is interesting, for although it is con- 
cerned with two-dimensional applications, it compares critically the 
performance of both finite difference and finite element algorithms. The 
computer time required by the finite element method is shovm to be much 
larger than chat required for a finite difference solution and this is 
identified to be due to the finite element processes of nume^'ical quadra- 
ture, assembly and solution. As has already been mentioned, the finite 
element practitioners who are currently working with three-dimensional con- 
figurations are aware of these restrictions and are developing techniques 
designed at reducing the CPU time penalty of the finite element method. 

Other criticisms of the finite element method, such as those concerned with 
the ability to include upwinding, have been removed by recent finite element 
developments described in the preceding sections. Cooke does emphasize that 
accurate results appear to be the rule rather than the exception with the 
finite element method, and that the variable grid capability is the method's 
greatest asset. 
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RECENT TRENDS IN SOLVING FINITE ELEMENT 
SYSTEM EQUATIONS 

In addicion to the improvement i in finite element modelling «K>rk 
already deacribed, baaic algorithmic development atudiea are currently being 
carried out to improve the competit iveneaa of the finite element aiethod. 

Thia work iu aimed at providing economic aolution techniquea for the ayatem 
equationa. Notable here ia the work of Hughea et al. (ref. 25) who propoae 
an element-by-element aolution acheme which avoida the aaaembly proceaa by 
working solely with the element equationa. Thia approach appeara attractive 
but it needa to be more thoroughly teated. On an operationa count baaia 
(see Fig. 3) it should certainly prove advantageous for 3D problems involv- 
ing a single degree of freei'om per node, but will it prove as competitive 
for systems of interest where the element matrix is typically 40><40, and is 
it a vectorizable process? Sample results of applying the method to a 
simple problem are shown in Figure 4. 

An alternative approach ie suggested by the work of Park (ref. 26) in 
which the actual equation system ia taken to be of the form 

(M + At K) ^ (6) 


where n denotes the time level. Matrix inversion can be avoided if 
equation (6) Is replaced by 


_1 _1 n 

M(^ + AtM L)(^ + AtM ];)♦ -_f 


(7) 


where 


T 

K - L ♦ U, L - U 


( 8 ) 


and ^ and ^ are lower and upper triangular matrices, respectively. 
Equations (6) and (7) agree to 0(A t^ ) but it was found that, although 
stability could often be guaranteed, unacceptable errors are produced, even 
for intermediate values of At. Attempts at using a non-symmetr ic splitting 
(i.e. lT a u) in equation (8) were more successful, but the production of 
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a general nechod for performing thia aplitCtng was not forthcoming. Park, 
therefore, retained the aynnuetric aplit and replaced equation (6) by 


(m ♦ At^m ♦ At(K ♦ Atjc)] ■ f^ (9) 

where m and are diagonal matricea yet to be determined. The equation 

(M ♦ At^m) { a ♦ At L) a ♦ At U)} ♦ ” ■ f (10) 


ia uaed inatead of equation (7) where now 


L - (M ♦ At^m)' (L ♦ — k) 
~ 2 


A 

u 


(M ♦ At^m) 



( 11 ) 


and it followa that equation (10) reduces to equation (6) provided that 

{m + k + (L ♦ k) (M ♦ At^m)’ (U ♦ — k)} - 0. (12) 

2 "*“ ~ ~ 2 

It is not possible to apply equation (12) at level n as required since 
is unknovm and, therefore, Jc and m are determined by applying 
equation (12) to the solution at level n-1. This approach has been shown 
to remove the poor accuracy associated with the first method for soaie simple 
structural dynamics examples, but no evidence is available %diich supports 
its adoption for more complicated problems. 

In fluid flow modelling it is frequently essential Co use a fine mesh 
in the vicinity of solid boundaries, to ensure adequate representation of 
the boundary layer, %«hereas a much coarser mesh may be used in the flow 
interior. The stability criterion associated wj*;‘» the fine mesh miy rule 
out the use of an explicit method, whereas the explicit stability limit for 
the coarse mesh alone could be acceptabl.*. The solution in this case may be 
the use of a mixed implicit-explicit technique in which the elesients forming 
the fine mesh are treated implicitly while those in the coarse mesh are 
created explicitly (ref. 27). As an example. Figure 5 shows a finite ele- 
ment mesh with 3 explicit elements end 2 implicit elements and the structure 
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of the resulting matrix which must be inverted assuming 2 unknowns per node. 
In the solution process the zero entries outside the profile need not be 
stored or operated upon. A finite difference equivalent of this type of 
approach has been investigated by Shang (ref. 28). He found an order of 
magnicude Increase in speed over the explicit method could be achieved, but 
numerical oscillations were observed with eome coarse grid configurations in 
3D. 


THE DEVELOPMENT OF A FINITE ELEMENT COMPUTER PROGRAM FOR 
ANALYSIS OF AERODYNAMIC HEATING IN THREE-DIMENSIONAL 
HIGH SPEED COMPRESSIBLE VISCOUS FLOW 


It is clear at the outset that any computer program for the local 
analysis of three-dimensional flows will require sophisticated pre- and 
pr,8t-processing software for displaying the initial configuration, boundary 
conditions and the computed results in a convenient fashion. It is not 
proposed, at least initially, to spend a large amount of time in developing 
such software but to utilize and develop existing and proposed facilities 
available at NASA/Langley, e.g. the pre-processing might be accomplished via 
a variant of PATRAN-G while the post-processing can utilize the software 
currently in use with, and under development for, 3D finite difference 
codes. The main effort can then be directed at producing a finite element 
computer program for solving the equations of high speed compressible flow 
in, eventually, three dimensions. As has been shown by the work of Gresho 
(refs. 17-19), the program will only be capable of handling the size of 
problem envisaged (of the order of 10® unknowns), provided it is highly 
vectorized so as to make full use of the vector processing speed of the 
Langley Cyber 203 (refs. 29-31). The optimimi finite element solution algo- 
rithm is by no means apparent at this stage and parallel algorithm and tech- 
nique investigation will be required before large scale problems can be 
attempted. The stages of development envisaged are shown in Figure 6. The 
2D code, which will be produced initially, will be written quite generally 
as far as the number of dimensions is concerned, so that the transfer to 3D 
at the end of Stage 3 should prove straightforward. It is envisaged that an 
essential feature of the development will be the comparison of numerical 
results with experimental observations. Preliminary discussions have taken 
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place with regard Co possible experiments chat might be useful for code 
validation in the initial stages of the development. In the later stages, 
code validation will become more critical and experimental support essen- 
tial. These later experiments would have to be detailed during Stage 3, by 
which time any uncertainties and manerical problems needing further investi- 
gation would have become apparent. 

CONCLUDING REMARKS 


This report reviews the current status of viscous flow modelling by the 
finite element technique. It is apparent that very little published finite 
element work is available in Che area of high speed compressible viscous 
flow and that this is a challenging field for the finite element method. A 
review has also been made of some of the techniques which are currently 
under investigation for the solution of finite element system equations. A 
proposed program of work leading to the production of a finite element simu- 
lator for local analysis of 3D high speed viscous compressible flow, 
governed by the full Navier-Stokes equations, has been outlined. Experimen- 
tal support has been identified as an essential feature of the simulator 
validation process. 
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APPENDIX 


Using tensor ni 
compressible viscous 

Continuity 

Momentum 

Energy 

Equation of state 
(perfect gas) 

Stress/rate of 
strain law 


>tation and the summation convention, the equations of 
flow can be written as 


do 9 

II + (pu.) - 0 

3t 3x. J 
3 


3 3 

— (pu .) + (u . pu . + p 6 . . - o. .) ■ 0 

3t ^ 3x. J ^ 

J 


i-J 


3 \ 3 / 1. 3T . 

— (pe) + (u. pe + u. p - or . u - k ) 

3t 3x. ^ ^ ^ 3x. 

J J 


e ■ C T p • pRT 

V 


o . . - M (T) . 
ij 


3u , ?u . 
J J- 


«(T) ®“k , 

ij 
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Figure 3. A comparison of operation counts between the element-by- 
element algorithm of Hughes et al. (ref. 25), and the 
standard globally implicit methods for a finite element 
solution of a problem involving N nodes with a single 
degree of freedom per node. 
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Figure 4. 


Element-by-e lenient algorithm used to solve the equation 
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wi th time step At = 0.025. 
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Figure 5. Two-d imens iona 1 impl i c i t -expl i c i t finite-element mesh and 
profile of system matrix (from Hughes and Liu, ref. 27). 
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Stage 1 

Aug. 82 - June 83 

The production of a 2D unsteady 
code VISC0H2D designed to be as 
simple as possible and initially 
to be tested on 10 shock tube 
problems. 


Evaluation of solution algorithms 
and identifying possible artificial 
viscosity models. 


Stage 2 

June 83 ~ Sept. 83 

2D calculations using VISC0K2D 
- numerical comparisons. 


Stage 3 

Sept. 83 ~ June 84 

Implementation and testing of 
Improved solution algorithms 
into VtSC0K2D - numerical and 
experimental comparisons. 


Stage 4 
June 84 - 


Evaluation of methods designed to 
improve code efflcie'icy t.g. 
integration techniques. 


Testing and running of the 
optimum 3D code VISC0H3D - 
numfTical and experimental 
compar i sons . 


Figure 6. 


Proposed stages in the development of 
computer code for solving problems of 
speed viscous cocg>ress ib le flows. 
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